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Frequency correlations in reflection from random media 
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We present a theoretical study of frequency correlations of light backscattered from a random 
scattering medium. This statistical quantity provides insight into the dynamics of multiple scattering 
processes accessible both, in theoretical and experimental investigations. For frequency correlations 
between field amplitudes, we derive a simple expression in terms of the path length distribution 
of the underlying backscattering processes. In a second step, we apply this relation to describe 
frequency correlations between intensities in the regime of weak disorder. Since, with increasing 
disorder strength, an unexplained breakdown of the angular structure of the frequency correlation 
function has recently been reported in experimental studies, we explore extensions of our model to 
the regime of stronger disorder. In particular, we show that closed scattering trajectories tend to 
suppress the angular dependence of the frequency correlation function. 

PACS numbers: (030.5620) Radiative transfer; (290.4210) Multiple scattering; (290.1350) Backscattering; 
(030.6140) Speckle. 


I. INTRODUCTION 

Correlations in frequency space have proven to be an 
important tool for analyzing the dynamics of multiple 
scattering processes of waves in disordered media. For 
fixed scattering particles, i.e., in a setup where the fre¬ 
quency dependence is not influenced by eventual move¬ 
ments of the scatterers [T], the second order frequency 
correlation function of multiply scattered field ampli¬ 
tudes can be represented as the Fourier transform of the 
time-of-flight distribution [2]. This relation has been ex¬ 
ploited, e.g., in early experimental investigations 
studying intensity-intensity correlations of light or mi¬ 
crowaves transmitted through random samples in depen¬ 
dence of the system parameters. In the same spirit, more 
recent experimental studies of intensity-intensity corre¬ 
lations in transmission demonstrated how these can be 
used to extract dynamical transport parameters like de¬ 
lay times EE! or the diffusion constant p3j in transmis¬ 
sion experiments through disordered media. 

This picture of the immediate connection between fre¬ 
quency correlations and the path lengths of the scattering 
sequences becomes even richer when considering the case 
of reflection from a random medium where, due to differ¬ 
ent contributing scattering processes as compared to the 
transmission scenario, the path length distribution de¬ 
pends crucially on the backscattering angle d. This fact 
manifests itself in the effect of coherent backscattering 
EHI23, leading to a sharp angular cone of the backscat¬ 
tered intensity around the direction $ = 0 corresponding 
to exact backscattering. Consequently, this dependency 
of the path length distribution on the backscattering an¬ 
gle leads to a characteristic angular peak of the width 
of the frequency correlation function M■ In a previ¬ 
ous work on reflection from two- and three-dimensional 
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random wave guides m , the narrow angular peak that 
results in the diffusive regime of weak disorder could not 
be resolved, due to the discrete scattering channel geom¬ 
etry. On the other hand, however, this work explicitly 
reports on effects of the transition between the localized 
and the diffusive regime on the delay time statistics. In 
particular, a significantly broader distribution of delay 
times is observed in the localized regime (i.e., for strong 
disorder) than in the diffusive regime (i.e., for weak dis¬ 
order). In a related discussion concerning the physics of 
the localization transition, correlation properties of radi¬ 
ation reflected from random media, both, in angular [16] 
or frequency m space, recently have been observed to 
undergo drastic changes as the disorder becomes suffi¬ 
ciently strong. Being directly related to the dynamical 
features of the underlying scattering processes, frequency 
correlations are considered promising for obtaining a bet¬ 
ter understanding of the change of scattering mechanisms 
and propagation properties at (or close to) the Ander¬ 
son localization transition. In that context, recent ex¬ 
perimental results m indicate that the characteristic 
angular peak of the width of the frequency correlation 
function mentioned above displays a drastic breakdown 
in the regime of stronger disorder near the localization 
threshold. The fact that this effect has so far not been 
explained theoretically provides the main motivation for 
the present paper. 

In this work, we focus on the case of reflection from 
a random scattering medium and present detailed the¬ 
oretical studies of the intensity-intensity correlations in 
the frequency domain. This represents in two ways a 
completion and an extension of the previous works de¬ 
scribed above: On the one hand, we enrich the knowl¬ 
edge about frequency correlations in reflection as com¬ 
pared to the case of transmission by tracing back the 
properties observed in experiment to the nature of the 
contributing scattering processes for a backscattering ge¬ 
ometry. In particular, we show that, for the case of a 
laser beam with large transverse width p £ (where i 
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denotes the mean free path) and for backscattering an¬ 
gles d > 0, frequency correlations are also affected by 
propagation outside the scattering medium. We derive 
an analytic expression for this dependency, which allows 
to map the frequency correlation function between field 
amplitudes onto the path length distribution of multi¬ 
ple scattering trajectories within the scattering medium 
(and vice versa). Furthermore, by investigating the im¬ 
pact of varying the amount of disorder in the system, 
we discuss how frequency correlations might help to shed 
more light on possible changes occurring at the localiza¬ 
tion transition on the level of the underlying scattering 
mechanisms. 

After building up the theoretical frame in Sec. [TTJ we 
show in Sec. m how to describe frequency correlations 
between multiply scattered field amplitudes in the regime 
of weak disorder, thereby developing a simple relation be¬ 
tween the frequency correlation function and the distri¬ 
bution of path lengths of the multi-scattering sequences 
contributing to reflection. In Sec. IV we apply this rela¬ 
tion to treat frequency correlations between field intensi¬ 
ties in the regime of weak disorder, providing a theoret¬ 
ical interpretation for previous experimental results |14] 
in this regime. Finally, in Sec.[Vj in order to investigate 
the behaviour of frequency correlations near the localiza¬ 
tion transition, we explore several corrections to the weak 
disorder approximations applied in Sec. m Whereas cor¬ 
rections to the Gaussian GA-approximation are shown to 
be negligible in the case of our backscattering setup, we 
demonstrate that recurrent scattering trajectories [TF] do 
indeed tend to suppress the angular peak of the frequency 
correlation function, as it has been observed in experi¬ 
ment Il4l. 


II. THEORETICAL FRAME 

In order to describe the frequency correlations of light 
reflected from a disordered medium in terms of the un¬ 
derlying multiple random scattering processes, we work 
in the frame of the following model: The disorder is de¬ 
scribed by a random potential V(r) = 6n(r)/(n) describ¬ 
ing relative fluctuations of the dielectric constant n with 
respect to its mean value ( n ) within the region of the scat¬ 
tering medium. The potential V (r) is modeled as a Gaus¬ 
sian random process, which is characterized by its mean 
value (V(r)) = 0 and correlation function (V(r)V(r')). 
Here, and in the following, angular brackets denote dis¬ 
order average over different realization of the disorder. In 
the regime where the wave length of the scattered wave A 
is much larger than the correlation length of the disorder 
potential, the correlation function can be approximated 
by: 

(U(r)U(r')) = u 2 5(r — r'), (1) 

where we introduced the pre-factor u 2 characterizing the 
scattering strength of the disorder. In scalar approxima¬ 
tion, a wave with wave number “fc = y propagat¬ 


ing in the presence of the disorder potential satisfies the 
scalar Helmholtz equation 


V 2 + “fc 2 (l + V(r)) 


“VO") = p{ r), 


( 2 ) 


where p( r) denotes a given source distribution. 

Instead of treating a special solution for one partic¬ 
ular realization of the disorder, we will be interested 
in average properties. First, we describe the average 
propagation between the points r and r' in space, i.e., 
we consider the average Green function LJ G( r — r') := 
f J Gy(r,r')) obtained by ensemble averaging the Green 
function “Gy(r,r') associated to Eq. ([2]) for one special 
configuration V of the disorder. The average propagator 
can be written as JBj 


1 r'| 

J G(r,r') = - —- jt— e~ 

4-7T r — r' 


(3) 


which exhibits, compared to the vacuum case V = 0, 
where the propagator is given by a spherical wave, ad¬ 
ditionally a term describing exponential damping. The 
damping constant t is referred to as the scattering mean 
free path , which describes the average distance between 
two successive scattering events. With t as the character¬ 
istic length scale of the disorder, we further introduce the 
disorder parameter kl classifying the amount of disorder 
in the system by comparing the wavelength of the scat¬ 
tered wave to the scattering mean free path. Hence, the 
disorder is said to be weak for large values kt 1 of the 
disorder parameter, whereas the strength of the disorder 
increases for decreasing values for kt. For weak disor¬ 
der, the mean free path turns out as t = 47r/(fc 4 u 2 ) [T8] . 
Moreover, in this regime, scattering of the average wave 
intensity can b e d escribed by ladder and crossed dia¬ 
grams (see Sec. Ill) which, in turn, are constructed from 
average Green functions, Eq. (0), and scattering events 
defined by the potential correlation function, Eq. 0 , in¬ 
troduced above. 


III. FREQUENCY CORRELATIONS BETWEEN 
FIELD AMPLITUDES 

We aim at describing the statistical properties of light 
of two different frequencies ui and w+Q which is backscat- 
tered from a random scattering medium. For the sample, 
we choose a slab geometry of thickness L, within which 
the disorder is characterized by the scattering mean free 
path l. The setup under consideration is sketched in 
Fig-0 In the following, we consider purely elastic scat¬ 
tering and we neglect the effects of absorption or internal 
reflections on the boundaries. 

The first quantity we focus on is the field correlator 

pm = r +n vww “poutm, (4) 

which can be seen as a generalized intensity composed of 
the two outgoing fields “ +n Vw($) and “ if out ($) which 
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FIG. 1. Scattering setup under consideration: Light of two 
different frequencies uj (orange) and uj + S2 (blue) is backscat- 
tered under the backscattering angle from a random scatter¬ 
ing medium where the disorder is characterized by the scat¬ 
tering mean free path t as the average distance between two 
successive scattering events of the multi-scattering sequence. 
Disorder average is taken over several different realizations of 
the random potential. 


FIG. 2. Scattering processes yielding minimal dephasing be¬ 
tween the wave w+n ipout at frequency u: + fl (solid lines) and 
its complex conjugate “V’out at frequency uj (dashed lines) 
as they are taken into account to compute the field correla¬ 
tion function see Eq. 0. within the approximation 

of weak disorder: Ladder (L) and crossed (C) process, where 
the same scatterers are visited in the same or in the reversed 
order, respectively. 


are scattered into the same direction parametrized by the 
scattering angle $, but differ in frequency by fl. To com¬ 
pute this quantity, we apply the following approximation, 
as it is commonly used in the treatment of disorder aver¬ 
aged intensity propagation in the case of weak disorder: 
We assume only those interference processes which lead 
to a minimum dephasing between the wave and its com¬ 
plex conjugate to yield non-negligible contributions on 
average. Minimal phase shift is achieved either if the 
wave and the complex conjugate visit the same scatter¬ 
ers in the same order (ladder propagation, L) or in re¬ 
versed order (crossed propagation, C). The equivalence 
of ladder and crossed processes in media exhibiting time- 
reversal symmetry for fl = 0 leads to the famous coherent 
backscattering effect png. Within this approximation 
of minimal dephasing , the ladder and crossed scattering 
processes shown in Fig. [5] have to be taken into account 
in order to compute F$(Q): 


using Eq. ([3]) in Fraunhofer approximation [l8j : 
1 _i“fe|R.—r| 

X-V / T-fc \ x ° -—5-- 


~ -,-tg(asin0-zcoBfl) 6 r — p 7 \ 

4 ttR ’ [ ’ 

where we used R = Re ou t and “k out = “fce OMt = 
^(cos$, 0, — sin#) to parametrize the outgoing scatter¬ 
ing direction. Hence, z/ cosd gives the distance travelled 
by the field inside the medium from the point of the last 
scattering event to the boundary surface, cf. Fig. [2j The 
wavy lines in Fig. [2] connecting the points iq and r 2 in¬ 
dicate the ladder (L) or crossed (C) propagation process 
of the average intensity between these two points, which 
we describe with the mixed frequency intensity propaga¬ 
tors u Pl{ ri,r 2 ) and n ?c(ri,r 2 ), respectively, where the 
superindex f l indicates a frequency shift f1 between the 
two propagating fields. Using the following notion of a 
single-step propagator 


F^n) = F^\n) + F^ c \ft), (5) 

where F^ L \Cl) and F^ C \Q) are constructed out of the 
following components: The incident field ,jj 'fi n (r) of am¬ 
plitude ifo is described by a Gaussian beam of width p 
and wave vector = “e z perpendicular to the x-y 

plane as the surface of incidence: 


n Po(r 1 ,r 2 ) 


w+n G(iq,r 2 ) aj G*(r 1 ,r 2 ) 

1 _ |ri- r 2 , 

(47r) 2 |iq - r 2 p e 


( 8 ) 


to describe a single scattering step between points iq and 
r 2 in space in terms of the average Green function given 
in Eq. ([3]) , the intensity propagator can be written in the 
following, self-consistent integral representation: 


“ifinir) = ip 0 e *e 2 *> 2 e 1 " 2 , (6) 

where r x = \J x 2 + y 2 denotes the transversal component 
of the vector r = (a :,y,z), i.e., its projection onto the 
plane of incidence. 

The outgoing radiation is described in the far field 
(R r, where R refers to the position of the detec¬ 
tor and r to a position within the scattering medium) 


n P L (ri,r 2 ) = u 2 S(r 1 -r 2 )+u 2 f dr n P 0 (r 1 , r)°P L (r, r 2 ), 

(9) 

where the term u 2 S(ji — r 2 ) represents single scatter¬ 
ing. Remembering the symmetry of the Green function 
under exchange of the spatial arguments, “G(rq,r 2 ) = 
u G(r 2 , iq), a single crossed scattering step, where one of 
the scattering paths is reversed compared to its ladder 
counterpart, is described by the same single-step prop¬ 
agator ^Pq (iq, r 2 ) as a ladder step. Hence, the crossed 
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sequence yields an equivalent propagation process - ex¬ 
cept for the fact that single scattering is excluded from 
the crossed propagator: 


n Pc(ri,r 2 ) = u 2 J dr n P 0 (ri,r) n P L (r,r 2 ). (10) 


Using the ingredients of Eqs. and ([9]), we identify 
the diagrams of Fig. [2] to be described by the following 
formulae: 


F^ L \n) = JJ dr ± dr 2 u+n ip in ( ri) u ^* n {r 1 ) n P L {r 1 ,r 2 ) 
x“+ n G out (r 2: nrG: ut (r 2 ,R), ( 11 ) 

= JJ dr lC lr 2 u+n ip in { r 2 ) “V4( r i) n Pc(ri,r 2 ) 

x^G^tfrr.RrG^^ra.R). (12) 


Inserting Eqs. ( |6|7| for the incoming and outgoing waves, 
we can explicitly write down the dependence on the fre¬ 
quency shift and the scattering angle t? as follows: 


F^ L \n) = 


i — R 
e „ n 


(rtl 2 


Pf } (U) = 


2 — pt 

e * 


dridr 2 e p 2 e A Zi + cos 2 t>) 

(13) 

_ »l+»2 (l 1 \ 

g 2 1 V v cosil ; 


(47rP) 2 

x n P L (n,r 2 ) pi,p 2 ,^) ) 

( r l~ ) 2 +( r 2~ ) 


dridr 2 e 


(4ttP) 2 

X e* "^((sci — X2)smi7+(« i — z 2 )(l—cosiT)) 

x n P c (r 1 ,r 2 )e i ? Ac(ri ’ r2 ’ ,J) , 


(14) 


where we introduced the enter/exit dephasing lengths 
(see below) for ladder and crossed diagrams: 


A L (r 1 ,r 2 ,^) = z 1 — x 2 sin d + z 2 cos t?, (15) 

Acr(ri,r 2 ,i?) = z 2 - x\ sini? + z\ cost?. (16) 


In Eqs. (p~3{ fl4|) , we recognize (apart from the constant 
term etwo terms depending on the frequency shift 
fl and thus giving rise to two different dephasing mecha¬ 
nisms: 

(i) Propagation within the scattering medium, (i.e., be¬ 

tween the first and last scattering event), described by 
the terms n PL,c( r u r 2 ). As evident from Eq. ([s]), a phase 
factor is picked up at each single step of the 

multiple scattering process. For the total sequence, these 
single phase factors add to a total phase e *T £ ( r i->r„) 
accumulated by the propagator n PL,c( r i,r n ), where 

—> r„) = I Y i ~ r i+il denotes the total path 

length of the scattering path from the initial scatterer at 
point ri to the final scatterer at r n (see also Appendix [A|. 

(ii) Propagation outside the scattering medium (i.e., 
before the first and after the last scattering event), de¬ 
scribed by the terms e l ~ XL ’ c ^ ri ' r2 ’^: These terms are 


defined by the positions ry 2 (in particular their x- and 
^-coordinates Xi >2 and zi i2 ) of the first and t he last scat¬ 
tering event, respectively, according to Eqs. PM - 


These considerations allow us 
in the following form: 


to rewrite Eqs. 


Ff’ C) (U)= JJ d£d\% ( 0 L ' C) (£,\)e i % ( - c+x+R \ (17) 


where 2l^ i,C ^(£, A) denotes the joint distribution of the 
path length £ and the enter/exit dephasing lengths A for 
waves scattered into direction t? in the case fl = 0 (see 
Appendix [AJ . The latter can be obtained by a numeri¬ 
cal Monte-Carlo simulation of a random walk within the 
scattering medium (see Appendix [Si) . 

Treating the dephasing lengths £ and A approximately 
as independent variables, which we may assume since 
A l,c are dominated by the transversal coordinate x, 
whereas the total pathlengtlr £ exhibits translational in¬ 
variance along the plane of incidence, allows us to factor¬ 
ize the joint distribution 

2l^’ C) {£, A) ~ a£’ C) (£) B ( ^ C) (A) (18) 


into separate distributions for £ and A, respectively. 

The integral over the dephasing length A of the en¬ 
ter/exit process can now be approximated analytically 
by the following steps: The transversal coordinates of the 
initial scatterer follow a Gaussian distribution of width 
p, whereas the longitudinal penetration depth is damped 
exponentially with damping factor £, cf. Eq. ([6]). Since 
we assume (which is in accordance with the exper¬ 

imental parameters in El and turns out to be a crucial 
condition to observe effects of crossed propagation pro¬ 
cesses as we exp ose in more detail in Sec. 0 , A l and Ac, 
see Eqs. ( fl5|l6l ), are hence dominated by the transversal 
coordinate x. We thus approximate Ac ~ Ac — x sin 
where x = r sin <p is expressed in polar coordinates, with 
uniformly distributed angle <p and Gaussian distributed 
radial coordinate, i.e., Pt.(r) = (y^p) -1 exp[— r 2 /p 2 ]. 
The integral over A in Eq. © thus gives rise to the 
following term: 

h(il)= [ dr P t (r) / ^e^sin^sintf 

Jo Jo 2 tt 


-i(Q)V sin 2 # T 
= e 8 K c ' p 1 n 


1(2) p 2 sin 2 $ 


(19) 


where Io(x) denotes the zeroth modified (hyperbolic) 
Bessel function of first kind. Thereby, the total corre¬ 
lation function /^(U), see Eq. ([ 5 ]), can be expressed in 
terms of the path length distributions A# ’ {£) of lad¬ 

der and crossed propagation sequences as follows: 

F*(n)~h{Q) f d£ [A^\£) + Af\£)Y^ c+R \ 

( 20 ) 

Conversely, a measurement of the field correlation func¬ 
tion Ffl(£l) can be used to determine the distribution 
of path lengths A by inverse Fourier transformation of 
Eq. (p0|). Note that, in the case of a Gaussian beam 
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FIG. 3. We compare the normalized frequency correlation 
function C^(f2)/C^(0) in Gaussian approximation according 
to Eq. (221 evaluated in the approximation of minimal dephas¬ 
ing (using ladder and crossed diagrams as shown in Fig. [2| 
obtained numerically (solid lines) to the semi-analytic descrip¬ 
tion in terms of the path length distribution, see Eq. ( |20| ) 
(dashed lines). Good agreement is obtained for the differ¬ 
ent backscattering angles d = 0.1 (black), 0.3 (blue) and 0.5 
(red) under study. The system parameters in this plot are 
given by the disorder parameter k£ = 5, the width of the in¬ 
coming Gaussian beam p = 100 1 and the optical thickness of 
the scattering medium b = L /1 = 30. 


FIG. 4. Normalized width of the frequency correlation func¬ 
tion Af2(i?)/Af2i(i?) as a function of the backscattering angle 
d for different values for the disorder parameter kl = 15 (red, 
solid line), 10 (green, dashed line), and 5 (blue, dot-dashed 
line) in weak disorder approximation. The data describes pre¬ 
vious experimental results well in the regime of weak disorder 
1 but fails to describe the breakdown of the peak of the 
frequency broadening for stronger disorder (i.e., ki ~ 5) re¬ 
ported in experiment [14] . The inset shows the width A Q.l (d) 
in the case where crossed diagrams are excluded. Parameters: 
optical thickness b = L/£ = 30 and beam width p = 1001!. 


with large width p and backscattering angles d > 0, 
also dephasing due to propagation outside the scattering 
medium must be taken into account through the factor 
h{ fi) in Eq. ((20). 

Comparable experimental measurements of the path 
length distribution for transmission through a random 
scattering sample have been performed in )19l . As the au¬ 
thors point out, however, a measurement of the intensity- 


intensity correlation function (see Sec. IV) - which is in¬ 
sensitive to the phase of the field correlation function 
F{)(ii) does not contain the full statistical information 
needed to directly recover the path length distribution 
without imposing additional assumptions on the latter. 
It is possible, though, to obtain the required phase infor¬ 
mation using third order intensity correlations (20ll22| . 


IV. FREQUENCY CORRELATIONS BETWEEN 
INTENSITIES 


Since photodetectors are sensitive to the intensity of a 
light field, most experiments measure frequency correla¬ 
tions between intensities rather than between field ampli¬ 
tudes. In the following, we will consider the normalized 
quantity: 


c tf (n) = 


^out W) 

<“ +n IoutW) ( ^ lout m 


-i. 


( 21 ) 


where u I ou t{d) = \ w ipout(f})\ . In order to calculate 
Eq. (211, we apply the following two approximations 


valid for weakly disordered media: First, we apply the 
Gaussian approximation , assuming the scattered fields 
“VwO?) to follow Gaussian statistics. As a consequence, 
the four-field average in the numerator of Eq. (21) factor¬ 


izes according to Wick’s theorem and can be expressed 
in terms of two-field contractions: 


C'tf(Q) ~ Cj(S2) = 


F*(0) 


( 22 ) 


where F^(fl) denotes the correlation function between 
field amplitudes, see Eq. Q. Here, we assumed Q -C to 
and thus (“ +n / O ut0?)) - ThutiP)) = ^(0) in the de¬ 
nominator of Eq. pll). Second, we use, as in Sec. Ill 


the approximation of minimal dephasing in order to ex¬ 
press F$(Q) in terms of ladder and crossed diagrams, see 
Eqs. Q5|13|14p . 

We start by testing the validity of the approximate 
formula, Eq. (201, derived in Sec. Ill Fig. [3] confirms the 
very good agreement between Eq. (20) and the e xact sum 
of ladder and crossed propagators given by Eqs. ( |5|13|14P 
for three different backscattering angles d. Furthermore, 
Fig. pi] reveals that the frequency correlations C^(f2) sen¬ 
sitively depend on the backscattering angle d. 

To further investigate this behaviour, we examine the 
width A Cl(d) of the correlation function Cp(fl), defined 
by Ctf(Afl) = C&(0)/2 as a function of the backscat¬ 
tering angle d for three different values of the disorder 
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parameter k £, see Fig. [4j In order to distinguish the in¬ 
fluence of ladder and crossed diagrams, Fig. [4] shows the 
normalized width Afi(i9)/Afix,(i9), where AQb(i 9) is the 
width of the frequency correlation function if only lad¬ 
der diagrams are taken into account. Whereas the latter 
is independent of k£ and exhibits a slight monotonic de¬ 
crease as a function of ■& (see inset), a clear peak of the 
normalized width observed in Fig. [4] arises due to crossed 
propagation processes. According to Eq. (20), this be¬ 
haviour can be interpreted in terms of the path length 
distributions A^' C \C): shorter paths lead to a broader 
frequency correlation function. In exact backscatter- 
ing direction (£) = 0), the path lengths for ladder and 
crossed processes are identical - except for single scat¬ 
tering which, as mentioned above, contributes only to 
ladder processes. For this reason, the normalized width 
starts slightly below the value 1 at $ = 0. At larger 
backscattering angles, longer path lengths contribute less 
to coherent backscattering (see also Appendix 0 , and 
the width of the frequency correlation function corre¬ 
spondingly increases. At the same time, however, the 
total weight of crossed processes (integrated over all path 
lengths) decreases, such that the normalized width again 
approaches the value 1 for very large backscattering an¬ 
gles. In between, a maximum is found at $ ~ 2 l/ ■ which 
approximately corresponds to the angular width of the 
coherent backscattering cone [23]. The height of the peak 
slightly decreases for stronger disorder (i.e., lower values 
of k£). This can be traced back to the factor h(Sl) in 
Eq. (20) describing dephasing due to p ropagation out¬ 
side the scattering medium (see Sec. Ill), which becomes 
more relevant for larger backscattering angles. 

For weak disorder (i.e., kl ~ 15), the above results 
agree well with the experimental measurements [14j . For 
stronger disorder (k£ ~ 5), however, a drastic breakdown 
of the height of the peak was observed in the experiment 
[l4j . which is much more pronounced than the slight de¬ 
crease predicted by Fig. [4j We note that also the au¬ 
thors of }14] are able to reproduce their measurements in 
the regime of weak disorder by theoretical calculations 
comparable to ours. In contrast to [14], however, our ap¬ 
proach avoids additional approximations (such as, e.g., 
the diffusion approximation), and thus provides a more 
general and transparent theoretical interpretation. 


V. CORRECTIONS TO THE FREQUENCY 
CORRELATIONS FOR STRONGER DISORDER 

This discrepancy between the weak disorder descrip¬ 
tion of frequency correlations and the experimental ob¬ 
servations in stronger disordered media calls for a study 
of possible extensions of the weak disorder theory. In 
this section, we first discuss corrections to the Gaussian 
approximation, which we find to be negligible for the 
setup under study. In a second approach, we extend the 
approximation of minimal dephasing by taking into ac¬ 
count further scattering scenarios beyond the ladder and 



FIG. 5. Exemplary scattering scenario where the intensity 
propagators cross (black box H). The crossing is described 
by a Hikami-box [53] H = Ha + Hg + He- 


the crossed processes. 


A. Corrections to the Gaussian approximation 


The factorization of the four-field average in the course 
of the Gaussian approximation described in Sec. |IV| im¬ 
plies the underlying intensity propagators to follow in¬ 
dependent scattering sequences that do not share any 
common scattering points. Such crossings between the 
intensity propagators, for which an exemplary trajectory 
is sketched in Fig. [5j will induce non-Gaussian correla¬ 
tions among the fields. We denote with G™(Q) a contri¬ 
bution to the frequency correlation function as defined 
in Eq. (21) with exactly n— 1 points of crossing between 
the two mixed-frequency intensity propagators. We first 
consider the case n = 2 of a single crossing as follows: 
In the same spirit as Berkovits describing angular cor¬ 
relations in reflection [23] , we describe the crossing with 
a Hikami-box H = Ha + Hb + He [M] as depicted in 
Fig. [5] and the propagation of the intensity using diffusive 
propagators. In this frame, it is possible to evaluate the 
maximum of the C 2 -contribution at zero frequency shift 
(Q = 0) analytically [251 : 


r 2 (Q — nl — 1^2 3 cost? 

^ ^ 49 (, k£) 2 (S /£ 2 )’ 


(23) 


where S = np 2 denotes the surface of the Gaussian beam 


of width p. According to Eq. (|23|), C$ is a second order 


term in an expansion of the frequency correlation func¬ 
tion in 1 /(k£). The crucial dependence, however, stems 
from the behaviour C J oc £ 2 /S <x (£/p) 2 . The appear¬ 
ance of this factor can be explained as follows: if two 
photons enter the scattering medium at randomly chosen 
points within the transverse area S £ 2 of incidence, the 
probability of their crossing will be proportional to 1/S. 
The same reasoning applies to cases with more crossings 
(n > 2), which are of higher order in l/(k£), but scale 
similarly with 1/S: Once a single crossing has occurred, 
the probability of further crossings will be independent 
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FIG. 6 . Example of a recurrent scattering trajectory where 
the intensity propagates in a closed loop. We explore the con¬ 
sequences of these trajectories for the frequency correlation 
function as an extension of the weak disorder approximation 
of minimal dephasing, where only the contributions from lad¬ 
der and crossed processes (see Fig. [2| are taken into account. 



of S'. As already mentioned in Sec. m we assume p>f, 
since, in order to observe influences from crossed pro¬ 
cesses, for each scattering path a reversed counter path 
must exist. This is possible only if, both, the positions 
1*1 and r 2 of the first and of the last scattering event are 
covered by the laser beam, see also Eq. (291. The condi¬ 
tion p £, however, together with S = 7 rp 2 proves con¬ 
tributions containing crossing points, i.e., all terms C£ 
where n > 1, to be negligible as compared to C^(0) = 1, 
see Eq. (22), for the systems under study. Therefore, 
the Gaussian approximation ss remains a valid 
description of the frequency correlation function also in 
regime of stronger disorder. 


B. Influence of closed loops 

We now explore the consequences of additional scat¬ 
tering scenarios beyond the ladder and the crossed pro¬ 
cesses of Fig. [2] for the frequency correlation function. 
The nature of the underlying scattering processes when 
the strength of the disorder approaches the threshold 
of Anderson localization has attracted attention recently 
in the context of experimental results |16j : Changes in 
the behaviour of the angular correlations are reported as 
a consequence of so-called recurrent scattering trajecto¬ 
ries, i.e., closed loop scattering scenarios where the ini¬ 
tial and the final scatterer fall close to each other. Such 
a scattering process is sketched exemplarily in Fig. [ 6 ] 
The authors report a significant enhancement of the con¬ 
tributions from recurrent scattering trajectories to the 
backscattering signal as the value of the disorder param¬ 
eter approaches the localization threshold. 

In our Monte-Carlo simulations (see Appendix [b|) , we 
include a certain fraction of loop propagation processes 
into the scattering sequences in order to investigate their 
impact on the backscattering signal and the properties of 
the frequency correlation function. We account for the 
length of the scattering loop A to be distributed as 1/A 2 , 
as it was observed in experiment m, and the closed 
loops to exhibit angular properties like single scatter- 


FIG. 7. Composition of the mean backscattered intensity 
7 ( 1 ?) = 7 rsW + 7cms($) with different amounts of contri¬ 
butions 7 «s(i?) from Recurrent Scatterings (RS, blue) com¬ 
pared to contributions 7 CMs(d) from Conventional Multiple 
Scatterings (CMS, red) exhibiting the characteristic coher¬ 
ent backscattering peak around d ~ 0. We quantify the 
amount of RS contributions via the parameter T = 7 Rsi^ = 
0)/ r ycMs($ = 0) = 0.23 (solid lines) and 0.57 (dashed lines). 
Contributions from single scattering events have been ex¬ 
cluded. Remaining parameters as in Fig. [3] 


mg events. The amount of contributions from closed 
loops as compared to the contributions from conven¬ 
tional multiple scattering events (ladder and crossed pro¬ 
cesses) in the backscattering signal is quantified by the 
parameter T = 7 rs{$ = 0 )/ 7 cms{^ = 0 ) as the ratio 
of mean backscattered intensity stemming from recur¬ 
rent scattering events (RS) and from conventional mul¬ 
tiple scattering processes (CMS), respectively, in exact 
backscattering direction 7 ? = 0. This principle is demon¬ 
strated in Fig. [TJ where we show the split-up between 
7 rsW and 7 cmsW for different values of T. Here, 
the backscattered intensity u I 0 ut.{S) is expressed as di¬ 
mensionless quantity (‘bistatic coefficient ’ D3I) defined 
as 7 (0) = 47ri? 2 /(7rp 2 / 0 ) “ I 0 ut{9) with incident intensity 
lo = IV'ol 2 - 

After thereby having quantified the amount of re¬ 
current scattering, its impact on the frequency corre¬ 
lation function is shown in Fig. [ 8 j where we plot the 
normalized width of the frequency correlation function 
Af2(i?)/Af2(0) for different values of T. We observe the 
peak height in the frequency broadening to decrease as 
the ratio of RS contributions in the backscattering sig¬ 
nal increases. Since, as already discussed in Sec. [TVj the 
peak can be traced back to the influence of crossed di¬ 
agrams, this behaviour is consistent with the decreasing 
height of coherent backscattering cones observed in [7] 
Note that in the simulations behind Figs. [7] and [ 8 ] contri¬ 
butions from single scattering events are excluded from 
the backscattering signal to allow better comparability 
with experimental results HI] (see also Appendix |b|) . 

These results give first evidence that the inclusion of 











Backscattering angle $/rad 


FIG. 8. Normalized frequency broadening AO(i9)/Afor 
different values of the parameter T = 0 (black, solid line), 0.23 
(red, dashed line), 0.57 (green, dashed-dotted line), and 1.36 
(blue, dotted line) defining the amount of recurrent scatter¬ 
ing trajectories to the backscattering signal, see Fig. [7] The 
peak in the normalized frequency broadening significantly de¬ 
creases as the amount of recurrent scattering increases. The 
inset shows the width AOi,($) in the case where crossed dia¬ 
grams are excluded. Remaining parameters as in Fig. [3] 


scattering trajectories beyond the ladder and crossed di¬ 
agrams indeed leads to a significant breakdown of the 
peak in the width of the frequency correlation function as 
a function of the backscattering angle i9, in a similar way 
as it was observed experimentally M- In our approach, 
we account for the occurrence of recurrent scattering tra¬ 
jectories m on a phenomenological level, introducing 
them ‘by hand’ into our Monte-Carlo simulations. In fu¬ 
ture work, it will be interesting to verify or check these 
results in the light of a more rigorous, consistent micro¬ 
scopic treatment of multiple scattering in the regime of 
stronger disorder. 


VI. CONCLUSION AND OUTLOOK 

We theoretically investigated frequency correlations of 
light reflected from a random scattering medium. We 
established a direct relation between the frequency cor¬ 
relation function and the distribution of path lengths of 
multiple scattering sequences, which renders the path 
length distribution of the scattering sequences contribut¬ 
ing in reflection directly accessible to experiment. Fur¬ 
ther, we extended these considerations to a description 
of intensity-intensity correlations in the regime of weak 
disorder where our model agrees with previous, experi¬ 
mental results m- For the regime of stronger disorder, 
in which the same experiments report an unexpected and 
unexplained behaviour of the frequency correlation func¬ 
tion, we investigated several extensions to the weak dis¬ 
order model, in order to explain the discrepancy between 
experiments and theory in this regime. In particular, our 


results indicate that closed, recurrent scattering trajec¬ 
tories alter the properties of the backscattering signal in 
a way which points towards the experimentally observed 
behaviour. 

Therefore, future studies of the properties of the fre¬ 
quency correlation function in disordered media should 
focus on confirming the impact of recurrent scattering 
events and describing their influence in more detail. On 
the side of experimental investigations, it would be of 
great interest to combine the experiments presented in 
0 and 0, in order to obtain joint information about 
the width of the frequency correlation function and the 
amount of closed loop trajectories contributing at a cer¬ 
tain value of the disorder strength. From a theoretical 
point of view, the starting point could be two-fold: On 
the one hand, we lack a rigorous theoretical derivation 
predicting the fraction T of closed loop trajectories as 
a function of the disorder parameter k£. This might be 
achieved using self-consistent approaches to localization 
EZ10- On the other hand, the diagrammatic represen¬ 
tation of the recurrent scattering trajectories themselves 
is not yet fully clarified. In particular, taking into ac¬ 
count scattering sequences that return to a nearby scat- 
terer (instead of exactly to the same scatterer, as we have 
assumed in the present paper) might lead to a different 
angular behaviour of the underlying path length distri¬ 
bution. 

Starting from this point, the field is open for future 
investigations of the scattering dynamics close to the 
threshold of Anderson localization. The change of cor¬ 
relation properties in the angular m or frequency 0 
domain have already proven to be a fruitful approach to 
obtain a better understanding of the underlying physi¬ 
cal properties. Among these, frequency correlations ap¬ 
pear especially promising due to their direct connection 
to dynamical properties and time scales. Therefore, fu¬ 
ture work on frequency correlations in random media will 
hopefully help to gain deeper and more profound insights 
about light transport and localization phenomena. 
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Appendix A: Path length distributions 


Here, we provide the exact definitions of the path and 
dephasing length distributions used in Eqs. ( fl7| ) and (201. 
They characterize the propagation of average intensity 
for zero frequency shift (i.e., 0 = 0). We start by ex¬ 
panding the propagator Pl for 0 = 0 into a multiple 
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scattering series formally obtained by iteration of Eq. ©: 

OO 

Pl( y u v f ) = u 2 S(ri - 17 ) +M 4 P 0 (r ?; ,r / ) + Z u 2n 


n =3 


n— 1 


x dr 2 ... dr n _i P 0 (r J , r J+1 ) (24) 

J .7=1 

where ri = 7 and r„ = 17 in the second line. This 
allows us to write the propagator P+ as an integral over 
an (unnormalized) path length distribution: 


P L (ri,r f ) = J dCP£ (ri,r f ;C) 


(25) 


where 

p l{ y u y P p ) =u 2 5(r. l -r f )6(£) 


+u 4 P 0 (r i , r/)<5 (£ - (7 - 17 1 


n— 1 


+z u2n / rfr 2 • ■ • n p 0 (7,7+1) 


n=3 


3 =1 


n—1 


<<5 £ ~ Z i rfe - rfc + i 


(26) 


fe=i 


Inserting this into Eq. (13), the joint distribution of path 
and enter/exit dephasing lengths defined by Eq. ( [T7| ) for 
the case of ladder propagation is obtained as follows: 


21l(£, A) = 




dridr 2 e p 2 e i^ Zl+ »*«) 


(4ttP ) 2 

x Pf( r 1 ) r 2 ;>C)(5(A- A z ,(r 1 ,r 2 ,i?)). (27) 


The corresponding expression 2l c (£, A) for crossed pro¬ 
cesses is obtained in the same way from Eq. (14). Atten¬ 


tion must be paid, however, to the fact that single scat¬ 
tering processes are missing in P^ ( 7 , 17 ; C) as compared 
to the ladder case, i.e., p c( r i- 17 ; £) = P£(ri,rf,£) — 
u 2 S(r.i — r f)S(C). Integration over the dephasing length 
A yields the reduced path length distributions: 


MQ = 


1 


A C (£) = 


(47rP ) 2 
xPf ( JT,r 2 ;£), 

1 


dridr 2 e 


(•• 1 ) 




dridr 2 e 


p 2 e 1 


+ 3 L ) 2 + (rj-)2 

Zp 2 . 


(4t rP ) 2 

xg i 1 —x 2 ) sin^+(21 — 2 2 )(1 — cost?)) 

xPc (it, r 2 ;£). 


From the above expressions, we first observe that, for ex¬ 
act backscattering i? = 0 and large beam widths p —> 00 , 
the expressions for Al{C) and Ac(C) are identical (ex¬ 
cept for the single scattering contribution). For non-zero 
backscattering angles d > 0 , a non-vanishing complex 
phase factor is present in the crossed amplitude Aq(C). 
Upon disorder average, this leads to the suppression of 


long light paths for the crossed processes with rising an¬ 


gle mentioned in Sec. Ill In a similar way, crossed pro¬ 
cesses are also suppressed if the width p of the laser beam 




between the first and the last scattering event in Eq. (291, 


P-Ll 

U 


(28) 

-^(i+ssb) 


(29) 


which is typically comparable to the mean free path 


Appendix B: Monte-Carlo simulations 

In this appendix, we sketch the Monte-Carlo algorithm 
used for our simulations of the backscattered radiation. 
In the regime of weak disorder, where we apply the ap¬ 
proximation of minimal dephasing described in Sec. uni 
we insert the representation of the propagator as given 
in Eq. ( [24]) into Eqs. ( 13|14 ) to obtain multi-dimensional 
integrals / dri... dr n over the positions of scattering 
events. These can be solved numerically by virtue of an 
iterative procedure in which each single propagation step 
is simulated by drawing random numbers from appropri¬ 
ate distributions determining first the initial penetration 
depth and subsequently the angular orientation and the 
length of each respective step [29}. Averaging over dif¬ 
ferent realizations of the disorder is then achieved by re¬ 
peating the algorithm numerous times. The number of 
iterations N for the simulations amounted to N = 10 7 
for the data presented in Figs. [ 3 ] and [J] N = 10 6 in Fig. [7] 
and N = 10 5 in Fig. [8] At these numbers, the statisti¬ 
cal error on the computed quantities was assured to be 
smaller than visible on the scale shown. 

This scheme is extended to the regime of stronger dis¬ 
order by further taking into account recurrent scattering 
(RS) trajectories. Implementing RS paths into the scat¬ 
tering sequences in a consistent and flux conserving way 
is achieved as follows: First, we define a weighting fac¬ 
tor p as the probability that propagation occurs along a 
closed loop. Choosing different values for p corresponds 
to varying T in Figs. [7] and [HJ In accordance with ex¬ 
periment SSI, the length of the loop is distributed as 
Prs{ A) = 2f/A 2 for A > 21 (and 0 otherwise, i.e., the 
minimal length of a loop is assumed as 2£). Further¬ 
more, the two fields (solid and dashed lines in Fig. [6]) can 
complete the loop either in the same sense (clockwise or 
counter-clockwise) or in the opposite sense, giving rise 
to equally contributing ladder and crossed diagrams. In 
total, the recurrent scattering contribution thus results 
by modifying ladder and crossed propagators as follows: 
p l ( r *,r/;£) = Pg{Yi,vf,C) = %u 2 6(r t -r f )P RS (A) and 
inserting these expressions into Eqs. (27) and ( p~7] ) . At 
the same time, the conventional multiple scattering con¬ 
tributions are multiplied by the factor 1 — p. In order 
to ensure flux conservation EDI, the conventional ladder 
component is multiplied by an additional factor which 
we determine numerically by integrating the flux over all 
angles. 

For comparison with experiments on light scattering, 
we furthermore account for the polarization degree of 
freedom as follows: single scattering is filtered out in the 
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helicity conserving and linear non-conserving detection 
channel, whereas scattering processes of higher order are 
assumed to be equally split into the respective conserv¬ 
ing or non-conserving channel (for linearly or circularly 
polarized light). Likewise, crossed propagation processes 
fully contribute to the helicity conserving or linear con¬ 
serving channels, whereas they are filtered out in the re¬ 
spective orthogonal channels. Naturally, the condition of 


flux conservation must hold for the sum of two orthogo¬ 
nal channels. Under these premises, the results of Fig. [7] 
refer to the helicity conserving channel (where crossed 
processes contribute and single scattering is filtered out). 
Likewise, the data for A fl(0) in Fig. [8] also refers to the 
helicity conserving channel, whereas Af Il(6) (inset) is 
evaluated in the linear non-conserving channel (only lad¬ 
der processes, single scattering filtered out). 
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